In [2]:
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import pygrib as pb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import glob
#Get data path
datapath = glob.glob("./21072000/M-*grb2")
print(datapath)
In [10]:
"""
利用basemap套件繪製圖形，在read_grib範例中，可以得知變數的位置，這邊畫雨量跟10米風，初始場沒有雨量值
需要的步驟，
1.讀取grib資料，取出所需資料，包含2米溫度、10米u風、10米v風、經緯度
2.投影相關設定值(lon_0, lat_0, lat_1, lat_2)
lon_0的key:LoVInDegrees
lat_0的key:LaDInDegrees
lat_1的key:Latin1InDegrees
lat_2的key:Latin2InDegrees
3.繪製圖形
這邊會先設定雨量值的colorbar
"""
cwb_data=['None','#9BFFFF','#00CFFF','#0198FF','#0165FF','#309901','#32FF00','#F8FF00','#FFCB00','#FF9A00',\
'#FA0300','#CC0003','#A00000','#98009A','#C304CC','#F805F3','#FECBFF']
clevs = [0,1,2,6,10,15,20,30,40,50,70,90,110,130,150,200,300,400]
cmap = mcolors.ListedColormap(cwb_data,'preecipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)
cnt = 0 # 只是為了抓取1測的固定值用
for grb in range(len(datapath)-1):
if cnt == 0:
forget = pb.open(datapath[grb]).select()[:][0]
analDate = forget.analDate.strftime("%y%m%d%H")
lon_0, lat_0 = forget["LoVInDegrees"], forget["LaDInDegrees"]
lat_1, lat_2 = forget["Latin1InDegrees"], forget["Latin2InDegrees"]
lats, lons = forget.latlons()
proj = Basemap(projection="lcc",resolution='l',lat_0 = lat_0 ,lon_0 = lon_0 ,lat_1 = lat_1, \
lat_2 = lat_2, llcrnrlat = lats[0,0], llcrnrlon = lons[0,0], urcrnrlat = lats[-1,-1], urcrnrlon = lons[-1,-1])
cnt = 1
basicvar1 = pb.open(datapath[grb]).select()[:]
basicvar2 = pb.open(datapath[grb+1]).select()[:]
fcst = str(basicvar2[0].forecastTime)
valDate = basicvar2[0].validDate.strftime("%y%m%d%H")
accrain = basicvar2[61]["values"] - basicvar1[61]["values"]
u10m_data = basicvar2[66]["values"]
v10m_data = basicvar2[67]["values"]
fig = plt.figure(figsize=(16,12))
proj.drawcoastlines(linewidth=0.9, color='slategrey')
meridians = np.arange(60.,170.,5.)#設定繪製經度線的起點、終點及間隔
parallels = np.arange(-10.,60.,5.)#設定繪製緯度線的起點、終點及間隔
proj.drawmeridians(meridians, labels = [0,0,0,1],color='grey',linewidth=0.8)#繪製經度線 label = [左,右,上,下]
proj.drawparallels(parallels, labels = [1,0,0,0],color='grey',linewidth=0.8)#繪製緯度線
ctf = proj.contourf(lons, lats, accrain, clevs[1:len(clevs)-1], cmap=cmap, norm=norm, extend='both',latlon=True)
gap = 20 #沒有用gap的話，會很慘
proj.barbs(lons[::gap,::gap], lats[::gap,::gap], u10m_data[::gap,::gap], v10m_data[::gap,::gap] ,\
length=5,barbcolor="navy", latlon=True)
cbar = plt.colorbar(ctf,orientation='horizontal',pad=0.05,fraction=0.05)
cbar.set_label("$^\circ$C",size=14)
plt.title("Initial Time: {}".format(analDate), loc="left",fontsize=16)
plt.title("Forecast Time:{}".format(fcst),loc="center",fontsize=16)
plt.title("Valid Time:{}".format(valDate),loc="right",fontsize=16)
plt.suptitle("In-Fa Typhoon",y=0.85,fontsize=20)
plt.show()
plt.clf()
